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ABSTRACT 

We search the high-latitude portion of the COBEf\ Differential Microwave 
Radiometers (DMR) 4-year sky maps for evidence of a non-Gaussian temperature 
distribution in the cosmic microwave background. The genus, 3-point correlation 
function, and 2-point correlation function of temperature maxima and minima are all 
in excellent agreement with the hypothesis that the CMB anisotropy on angular scales 
of 7° or larger represents a random-phase Gaussian field. A likelihood comparison of 
the DMR sky maps to a set of random-phase non-Gaussian toy models selects the 
exact Gaussian model as most likely. Monte Carlo simulations show that the 2-point 
correlation of the peaks and valleys in the maps provides the greatest discrimination 
among the class of models tested. 

Subject headings: cosmic microwave background - cosmology: observations - 
methods: statistical 
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1 Introduction 



The angular distribution of the cosmic microwave background (CMB) reflects the 
density distribution at the epoch of last scattering. The statistical distribution of 
CMB temperature fluctuations is thus an important probe of the initial conditions 
for structure formation. Causally-connected regions at the surface of last scattering, 
as viewed from the present epoch, subtend an angle ~ 2°. Most of the current 
models for structure formation predict the CMB temperature fluctuations to be 
nearly Gaussian on larger angular scales. In inflationary cosmologies, this Gaussian 
distribution reflects the quantum origin of the density perturbations, while in 
topological defect models, it results from the superposition of a large number of 
non-Gaussian perturbations at smaller angular scales (the central limit theorem). 

The CMB temperature field may be conveniently represented in terms of 
spherical harmonics, 

T(M)=E a ^™(M) (1) 

£m 

and is said to be Gaussian if the distribution each coefficient a^ m follows a Gaussian 
probability distribution, 

PM = -j^= exp(-^|) 

with phases uniformly distributed in the range [0, 2n]. Gaussian fields have the 
desirable property that they are completely specified by the power spectrum or its 
Fourier transform, the 2-point correlation function. Thus, if the CMB is demonstrated 
to be a Gaussian field, not only have current models passed an important consistency 
test, but we are also spared the task of tabulating the full set of higher-order 
correlation functions. Deviations from Gaussian statistics would necessitate a 
dramatic reappraisal of current models. 

The COBE DMR experiment mapped the CMB at 7° angular resolution. A 
number of authors have tested the DMR one- and two-year sky maps and found good 
agreement with a Gaussian distribution (Smoot et al. 1994; Luo 1994; Hinshaw et 
al. 1994, 1995; Torres et al. 1995). Kogut et al. 1995 compared the two-year DMR 
maps to non-Gaussian toy models and found that the exact Gaussian model was the 
most likely descriptor of the data. In this Letter we compare the final 4- year DMR 
maps (Bennett et al. 1996, Kogut et al. 1996) to both Gaussian and non-Gaussian 
toy models and discuss the relative ability of various tests to discriminate between 
the Gaussian and non-Gaussian models tested. 
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2 Comparison with Gaussian models 



Cosmic variance, sample variance, and instrument noise can combine to create non- 
Gaussian realizations from initially Gaussian parent populations, while the central 
limit theorem pushes the CMB distribution back to Gaussian on large angular 
scales (Scherrer & Schaefer 1995). We infer the statistical properties of the CMB 
parent population using Monte Carlo techniques, in which a certain statistic (3-point 
correlation function, extrema correlation function, or genus) is calculated for many 
simulations of noisy CMB skies. The simulations include all relevant aspects of the 
CMB model, instrument noise, and data handling, including removal of the Galactic 
plane and fitted dipole/quadrupole subtraction. The same statistic is then calculated 
for the DMR data and compared to the range observed in the simulations. If the 
DMR data fall comfortably within the range of equivalent Monte Carlo realizations, 
we accept the model as in agreement with the data; otherwise, we reject the model 
at some specified confidence level. 

Unless otherwise noted, we model the CMB as a scale-invariant Gaussian field 
{P{k) oc Qfms-PS k n w hh n — 1) with quadrupole power-spectrum normalization 
Qrms-ps = 18 /iK (Gorski et al. 1996), evaluated at multipoles 2 < £ < 40 
filtered through the DMR window function (Wright et al. 1994). We compare 
these simulations to the DMR 4- year 53 GHz (A+B)/2 sum map. All statistics 
are evaluated only for the high-latitude sky (3881 pixels at |6|>20° and custom 
cutouts near Ophiuchus and Orion, Bennett et al. 1996) from which a fitted dipole 
and quadrupole have been removed. 

2.1 3-Point Correlation Function 

A Gaussian field is fully characterized by its two-point correlation function. All odd 
moments are identically zero, while all higher even moments may be specified in terms 
of the 2-point function. The 3-point correlation function and its zero-lag value, the 
skewness, are a first test for non-Gaussian signatures in the DMR data. 

The general 3-point correlation function is the average product of three 
temperatures with a fixed relative orientation on the sky: 

C 3 (e 1 ,9 2 ,0 3 ) = (T(hi)T(n 2 )T(h 3 )), 

where hi ■ h 2 = cos^i, h 2 ■ h 3 = cos9 2 , and h 3 ■ hi = cos# 3 . For computational 
reasons we evaluate only two special configurations, the most sensitive of which is 
the pseudo-collapsed case in which hi is nearly parallel to h 2 . Although the exact 
collapsed case (hi = h 2 ) is the simplest configuration, it involves terms quadratic in 
the observed temperatures and thus suffers from a noise bias (Hinshaw et al. 1994). 
We avoid this problem by using the map geometry to our advantage. Since the sky 
map pixels are 2?6 across and the DMR beam has a 7° FWHM, the nearest neighbor 
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pixels will have correlated sky signal but uncorrelated instrument noise. We thus 
define the pseudo-collapsed 3-point correlation function as 

C3 PC) = W i W j W k TiTjTk/ £ WiWjWk 
i,j,k i,j,k 

where Tj is the temperature in pixel i, w is the pixel weight, the sum on j runs over 
all pixels that are nearest neighbors to i, and the sum on k is over all pixels (except j) 
within a fixed angular separation of i. We also evaluate the equilateral configuration, 
Cf } , for which 0x = 6 2 = 9 3 . Since most pixels have 8 nearest neighbors, the pseudo- 
collapsed configuration sums over 8 times as many pixel triples as the equilateral 
configuration, resulting in a significant increase in sensitivity. 

Figure [TJ shows the pseudo-collapsed and equilateral 3-point correlation 
functions for the 4- year DMR 53 GHz (A+B)/2 maps, evaluated with unit pixel 
weight (w = 1) in 2?6 angular bins. The gray band in each plot shows the 68% 
confidence region derived from 2000 random realizations of scale-invariant CMB 
anisotropy superposed with instrument noise, and includes the effects of cosmic 
variance. Deviations from zero in C3 for either configuration are well within the 
range expected for a single realization of a stochastic process. We quantify this using 
a x 2 statistic, 

X 2 = E (A, " (3 a )) (D P - (S p )) (2) 

where D a is the DMR 3-point correlation function, (S a ) is the mean correlation 
function from the simulations (in this case identically zero), and 

M aP = 1 £ (Sg> - (S a )) (5f - (S fi )) (3) 

i 

is the mean covariance matrix between the angular bins a and f3 computed from the 
N = 2000 correlation functions in the simulations. We calculate 2000 such \ 2 
values by replacing the 3-point correlation function from the DMR data, D a , with 
the correlation function from the zth simulation, . The DMR maps fall near the 
median of the resulting distribution: 66% of the simulations had larger \ 2 than the 
DMR pseudo-collapsed 3-point correlation function, while 31% of the simulations had 
larger y 2 than the DMR equilateral configuration. The 3-point correlation function 
shows no evidence for a non-Gaussian CMB distribution. 

2.2 Genus 

Topology provides another test of the Gaussian hypothesis. The statistical properties 
of the CMB can be characterized by the excursion regions enclosed by isotemperature 
contours. The genus is the total curvature of the contours at fixed temperature, and 
may loosely be defined as the number of isolated high temperature regions (hot spots) 
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minus the number of isolated low temperature regions (cold spots). For a random 
Gaussian field, the genus per unit area is 



fa) = j^ei - e 

(Gott et al. 1990), where v is the temperature contour threshold in units of the 
standard deviation a of the field, and 8 C is the correlation angle, related to the 2- 
point correlation function 6*2(6') by 



C 2 (0) d9 2 |W - U ' 

The genus per unit area is a locally invariant quantity and is insensitive to incomplete 
sky coverage (e.g. removal of the Galactic plane). 

The genus of a Gaussian field depends on two parameters: the correlation 
angle 9 C and the standard deviation a. The DMR noise is nearly uncorrelated 
(Lineweaver et al. 1994); its correlation angle is determined by the pixelization. The 
correlation angle of a scale-free CMB depends only on the index n. For a noisy map, 
however, both the correlation angle and the standard deviation depend on the relative 
amplitude of the signal and noise terms and hence acquire a dependence on the CMB 
amplitude normalization Q rms -ps- Since 6^0^ <C # c ,cmb, the sensitivity of the genus 
to the CMB terms can be enhanced by smoothing the map, provided the smoothing 
dispersion remains less than # Cj cmb- We evaluate the genus at 3 smoothing angles with 
full width at half maximum 0° (no smoothing), 5°, and 10°. At each smoothing angle, 
a nearest-neighbor algorithm counts the number of hot and cold spots (contiguous 
pixels above or below the temperature threshold) at 11 thresholds v from -2.5a to 
+2.5cr in steps of 0.5a. We adopt the difference in spot number as an estimate of the 
genus and use the Monte Carlo simulations to calibrate any difference induced by the 
Galactic cut between this definition and the total curvature. The two definitions are 
identical in the absence of a Galactic cut, and the nearest-neighbor algorithm allows 
significant computational savings. 

Figure |2] shows the genus of the 53 GHz (A+B)/2 map smoothed with a 
Gaussian of 5° full width at half maximum (FWHM). The DMR genus is well within 
the range found by the simulations. As with the 3-point correlation function, we 
quantify this using a \ 2 statistic (equations and |3| with vectors D a and S a now 
consisting of the genus in k = 33 bins for the 11 thresholds and 3 smoothing angles 
evaluated simultaneously). The DMR genus is near the median of the Gaussian 
models: 51% of the simulations had larger y 2 . 

Since the instrument noise properties are known, we may interpret the genus 
in terms of the CMB model parameters. We evaluate the \ 2 °f the DMR genus on a 
2-dimensional grid in Q rms -ps and n, use the value at each grid point to determine 
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the likelihood 



C(Q. 



rms— 



ps,n) 



(2k) 



fc/2 exp(-|x 2 ) 
^/det(M) ' 



and search for the maximum in the resulting likelihood distribution. We also perform 
a probabilistic search by converting the x 2 & t each grid point to a probability 
P(x 2 ), defined as the fraction of simulations with x 2 larger than the DMR value 
at that grid point. We calibrate both methods and search for biases by generating 
2000 simulations with fixed parameters Q rm s-ps and n, fitting each realization, and 
examining the "cloud" of fitted values. The x 2 probability test has a small bias of 
0.8 /iK toward smaller values of Qrms-ps] the maximum likelihood test is unbiased. 

Both methods yield consistent results on the DMR 53 GHz (A+B)/2 data. 
The likelihood shows a strong ridge along the line Q rm s-ps — (17.3 ±1.8) — (8.8 ± 
0.2) (n—1) /iK with likelihood maximum at Q rms -ps — 24.5±5.0 /iK and n = 0.2±0.6 
(68% confidence interval). The value at the maximum has P(x 2 )=77%, and is not 
statistically preferred over the more precise value derived from the power spectrum 
of the combined DMR maps (Gorski et al. 1996). 

2.3 Extrema Correlation Function 

Another test of the statistical distribution is the angular correlation of the extrema 
(peaks and valleys) in the temperature field, defined as those points for which VT = 0. 
For a pixelized map, this reduces to the collection of pixels hotter or colder than all 
of their nearest neighbors. Specifying pixels hotter than their neighbors produces a 
set of "peaks" , while specifying colder pixels produces "valleys" . 

The number density of peaks or valleys is dominated by the noise properties 
of the map (Kogut et al. 1994). The clustering of the extrema, as measured by the 2- 
point correlation function of the maxima and minima, provides additional information 
on the underlying CMB temperature field. We define the extrema correlation function 
as 



(Kogut et al. 1995), where the sum over pixel temperatures is restricted to the 
set of maxima and minima pixels. Bond & Efstathiou (1987) provide analytic 
approximations for C cxt (6>) for random Gaussian fields but do not explicitly include 
the effects of instrument noise. Since the correlation properties of the non-uniform 
noise in the DMR maps are different from the underlying CMB temperature field, 
we smooth the maps with a 7° FWHM Gaussian prior to collating the extrema as a 
compromise between suppressing noise and removing CMB power at small scales. We 
use Monte Carlo techniques as described above to calibrate the resulting distribution 
of simulated correlation functions. 



Cext(#) 



J2ij WjWjTjTj 
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The clustering of extrema depends on the threshold v. We evaluate C ex t(0) 
at thresholds \v\ = 0, 1, and 2. Since, by definition, two peaks can not be adjacent, 
we ignore both the bin at zero angular separation and the first non-zero bin in all 
subsequent analysis. Simulations show that the results are dominated by the first 
few remaining bins; consequently, we speed processing by truncating the correlation 
function at separation 9 = 60° for a total of 22 angular bins at each of the three 
thresholds. 

The angular distribution of peaks and valleys in the DMR maps is in excellent 
agreement with Gaussian statistics. We calculate a x 2 value for the DMR extrema 
correlation function and compare it to 2000 simulations of CMB models with 
Qrms-ps — 18 /iK and n = 1. 70% of the simulations had larger x 2 than the DMR 
maps. 



3 Comparison with non-Gaussian models 

Figure |3] summarizes the comparison of the 53 GHz (A+B)/2 four-year sky maps with 
Gaussian CMB models. The x 2 derived from the 3-point correlation function, genus, 
and extrema correlation function of the DMR data fall near the median of the x 2 
distributions derived from the simulations: the DMR data are in excellent agreement 
with the hypothesis that the CMB is a realization of a parent population of Gaussian 
density perturbations. 

Given the large DMR beam, the agreement of the DMR data with Gaussian 
models does not necessarily place tight limits on non-Gaussian models. There are 
an infinite number of non-Gaussian distributions, from which we wish to select 
those representative of physically interesting physical processes. Ideally, we would 
use Monte Carlo realizations of specific cosmological models such as texture or 
global monopoles. Unfortunately, no analytic expressions exist for the temperature 
anisotropy produced by these models; generating a single realization requires a 
substantial computational investment to follow the causal field ordering in some large 
volume. Instead, we will modify the standard Gaussian model to test generic but 
interesting non-Gaussian amplitude distributions. 

We use a likelihood analysis of the statistical tests described above to compare 
the DMR maps against a set of toy models for which the spherical harmonic 
amplitudes ai m (Eq. [I]) are random variables with zero mean drawn from x% parent 
populations with N — 1, 5, 15, 30, or 60 degrees of freedom. We fix the variance cxf 
of the parent population at the same value as the standard Gaussian model, 

2 _ 4vr r[/ + (n-l)/2] r[(9-n)/2] 

1 ~ Wrms ~ ps) 5 r[Z + (5-n)/2]r[(3 + n)/2] 

(Bond & Efstathiou 1987), providing a simple, computationally efficient modification 
to the standard Gaussian model. Although the non-Gaussian amplitude distributions 
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tested here are skew-positive, the resulting temperature distributions (sky maps) are 
the convolution of the a,£ m with the spherical harmonics Y^ m and are thus characterized 
not by the skewness but by a positive kurtosis in the distribution of temperatures 
T (e.g., higher "wings" than a Gaussian distribution). The toy models provide a 
smooth transition from strongly non-Gaussian (N=l) to nearly Gaussian (iV=60), 
and may be compared to specific models of structure formation (Weinberg & Cole 
1992, Kogut et al. 1995). Cosmological models with rare high-amplitude peaks, 
typified by topological defect models such as texture or global monopoles can be 
compared to the x 2 models with few degrees of freedom, which tend to produce such 
features. 

Figure |] shows the likelihood of the 4-year DMR maps for the Gaussian and 
X% models, derived from the extrema correlation function. The likelihood is greatest 
for the exact Gaussian model, and decreases as the toy models become increasingly 
non-Gaussian. The Gaussian model remains the most likely descriptor of the data 
when the genus replaces the extrema correlation in the likelihood analysis. 

Given that the DMR data "prefer" the Gaussian model, how confident are we 
that the CMB is not well described by one of the \ 2 toy models? Since cosmic variance 
can produce a nearly Gaussian sky from a non-Gaussian parent population and 
vice-versa, marginalizing over model parameters can produce misleading confidence 
intervals (see the discussion of this point in Kogut et al. 1995). We quantify the ability 
of each statistical test to discriminate between Gaussian and non-Gaussian models 
using Monte Carlo simulations. We take 2000 simulations from each of the models 
in turn, and for each simulation generate the likelihood function £(Q rms _ps, N). We 
then count the number of times the Gaussian model is selected as "most likely." We 
obtain the DMR result (Gaussian model chosen as most likely) five times more often 
when the input model is Gaussian than when the input is any of the non-Gaussian 
models. The method is unbiased: if we ask how often the Xi5 model would be 
selected as most likely, we find the greatest probability if the input model is Xi5; and 
so forth. In this posterior sense, the CMB is five times more likely to be a realization 
of a Gaussian process than any of the \ 2 models tested. Table |I] gives the relative 
probability to observe a Gaussian sky if the CMB is a sample drawn from each of the 
toy models, derived using either the genus or the extrema correlation function (since 
all of the models tested have zero skewness, the 3-point correlation function does not 
discriminate among them). The extrema correlation function discriminates among 
the toy models twice as well as the genus. 



4 Conclusions 

We have tested the statistical distribution of temperature fluctuations in the DMR 
4- year 53 GHz (A+B)/2 summed sky map against both Gaussian and non-Gaussian 
models of CMB anisotropy using three statistical tests: the 3-point correlation 
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function, genus, and 2-point correlation function of temperature maxima and minima. 
A goodness-of-fit test using the x 2 °f each function shows the DMR data to lie near 
the median of the distribution of \ 2 values evaluated for randomly generated Gaussian 
models with amplitude normalization Q rms -ps = 18 //K and power-law index n — 1. 
The DMR data are thus in excellent agreement with the hypothesis that the observed 
CMB anisotropy represents a realization of a Gaussian parent process. We test 
an alternate hypothesis, that the CMB represents a realization of a non-Gaussian 
process, by evaluating the likelihood of the DMR data against a grid of non-Gaussian 
toy models whose spherical harmonic coefficients ae m are drawn from Xn distributions 
with N — 1, 5, 15, 30, or 60 degrees of freedom. The likelihood reaches a maximum for 
the exact Gaussian model. Cosmic variance produces significant cross-talk among the 
class of models tested and precludes the simple identification of confidence intervals 
on N. Simulations in which a similar likelihood analysis is performed on random 
realizations of either Gaussian or non-Gaussian skies demonstrates that the method 
correctly identifies the input model with a ~15% error rate against any of the other 
models; the error rate is nearly independent of N. The DMR result (Gaussian model 
selected as most likely) occurs five times more often when the input to the simulations 
is, in fact, Gaussian, than when it is any of the non-Gaussian toy models. We thus 
conclude that, not only do Gaussian power- law models provide adequate description 
of the CMB anisotropy on large angular scales, but that the non-Gaussian models 
tested are five times less likely to describe the true statistical distribution. 



We gratefully acknowledge the dedicated efforts of the many people responsible 
for the COBE DMR data: the NASA office of Space Sciences, the COBE flight team, 
and all of those who helped process and analyze the data. 
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Table 1: Relative Probability To Select Gaussian Model" 



Input 


Statistical Test 




Model 


Extrema Correlation 


Genus 


Gaussian 


5.1 


2.7 


xlo 


1.1 


1.0 


X30 


1.0 


1.2 


xi 5 


1.2 


1.2 


xl 


1.1 


1.3 


xl 


1.1 


1.1 



a Number of realizations from each input model for which the Gaussian model was 
selected as most likely, divided by the smallest such value to yield the relative 
probability of obtaining the DMR result from each input model (see text). 
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Figure 1: Three-point correlation function of the 4-year 53 GHz (A+B)/2 sum 
map, compared to simulated maps with Q rms -ps = 18 /iK, n — 1, and 53 GHz 
noise properties (thermodynamic temperature units). The error bars represent the 
uncertainty due to instrument noise alone, while the gray bands represent the 68% 
confidence region from the combined noise and Gaussian sky signal, (top) Equilateral 
triangle configuration, (bottom) Pseudo-collapsed configuration. The fluctuations 
about zero correlation are too large to be ascribed to instrument noise alone, but are 
consistent with the range of fluctuations expected from a Gaussian CMB. 
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Figure 2: Genus curve for the 4- year 53 GHz (A+B)/2 sum map, compared to 
simulated maps with Q rm s-ps — 18 A^Kj n = 1, and 53 GHz noise properties. The 
gray band represents the 68% confidence region from the simulations. The genus is 
in excellent agreement with the Gaussian CMB model. 
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Figure 3: x 2 °f t ne DMR maps compared to the distribution derived from simulations 
with the Gaussian CMB model, divided by the number of bins used to evaluated the 
X 2 ■ (top) Three-point correlation function, (middle) Genus, (bottom) Extrema 
correlation function. The x 2 values from the DMR data (vertical bars) are near the 
median of the distribution expected from Gaussian CMB anisotropy. 
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Figure 4: Likelihood function of the 4 year DMR 53 GHz (A+B)/2 extrema 
correlation function (marginalized over the CMB normalization Q rms -ps) f° r 
Gaussian and non-Gaussian models. 
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